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SECTION  1 
INTRODUCTION 


High-power  diodes  envisioned  for  the  1990s  pose  a  significant 
challenge  for  conventional  nunerical  simulation  techniques.  One  problem 
is  due  to  enhancement  of  magnetic  field  effects  in  the  central, 
high-energy  density  region  of  the  diode. 


in  the  conventional  simulation,  the  fundamental  electromagnetic 
frequency  (u>f)  dominates  the  choice  of  spatial  grid, 

€  2n c  ,  ( 1 ) 


to  achieve  spatial  resolution  of  the  wavelength.  Hence  the 
Courant-limited  electromagnetic  time  step  can  be  related  to  the  highest 
frequency  by 


ojf,6t  «  2n  . 


(2) 


in  most  regions  of  a  puised-power  system,  the  cyclotron  frequency, 


co  6t 
c 


eB6t 

m 


(3) 


presents  a  less  restrictive  constraint.  That  is,  a  conventional  particle 
kinematics  algorithm  can  usually  resolve  the  orbital  motion  using  a  time 
step  dictated  by  electromagnetic  wave  requirements.  However,  due  to  field 
compression  and  charge  leakage,  the  magnetic  field  in  the  central  region 
of  a  diode  can  be  sufficiently  large  that  orbita1  motion  is  not  resolved. 
Section  2  of  this  report  will  demonstrate  this  effect  upon  calculated 
particle  trajectories. 
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To  address  this  problem,  a  particle  subcycling  algorithm  was 
developed  and  implemented  in  MAGIC^.  The  basic  idea  of  subcycling  is  to 
integrate  the  electron  trajectories  using  a  time  step  smaller  than  the 
electromagnetic  time  step,  thus  permitting  trajectory  resolution  without 
an  excessive  nunber  of  electromagnetic  calculations.  To  test  variations, 
the  new  algorithm  in  MAGIC  contained  five  options  implemented  in  a  new 
executive  module.  These  options  and  the  executive  module  are  described  in 
Section  3. 


Tests  of  the  algorithm  were  performed  on  a  diode  simulation 
problem  based  on  Blackjack  5.  These  tests  included  a  detailed  analysis  of 
computational  time  broken  down  by  subroutines.  These  tests  demonstrated 
that  subcycling  did  yield  an  advantage  in  computational  speed.  However, 
tests  also  revealed  the  presence  of  numerical  instability  whenever  the 
speed  advantage  became  significant.  Upon  investigation,  we  were  able  to 
show  theoretically  and  in  simulations  that  subcycling  produces  instability 
under  the  conditions, 

u)  / oj  <  1  and  6tto  >  n/8  ,  (4) 

pc  c 

where  is  the  plasma  frequency  and  u)c  is  the  cyclotron  frequency. 

One  distinguishing  feature  of  this  instability,  termed  the  subcycling 
cyclotron  (SC)  instability,  is  a  nonphysical  growth  in  energy.  The 
subcycling  cyclotron  instability  is  discussed  in  Section  4  of  this  report. 

in  summary,  our  present  belief  is  that  subcycling  will  not  be  useful 
for  pulsed-power  applications.  However,  in  other  plasma  physics 
applications,  it  might  result  in  a  factor  of  two  or  more  computational 
efficiency.  To  achieve  such  results,  it  may  be  necessary  to  subcycle 
electrons,  while  retaining  a  larger  time  step  for  the  more  massive  ions. 


^B.  Goplen,  R.  Clark,  3.  McDonald,  G.  Warren,  and  R.  Worl,  "MAGIC  User's 
Manual,"  Mission  Research  Corporation  Report,  MRC/WDC-R-126,  April  1987. 
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SECTION  2 

SINGLE  PARTICLE  ORBITS 


2.1  PROBLEM  DESCRIPTION. 


The  problem  of  simulating  self-consistent  particle  motion  in  regions 
of  hiqh  magnetic  field  has  two  related  aspects.  First,  the  conventional 
kinematics  algorithms  will  not  correctly  compute  the  particle  orbit  if  the  time 
step  is  too  large.  Thus,  the  basically  circular  motion  inherent  in  physical 
drifts  (ExB,  BxVB,  etc.)  will  not  be  correctly  resolved.  Second,  charge-  and 
current-density  contributions  from  the  particle  orbit  will  not  be  correctly 
resolved  if  the  spatial  grid  is  too  large. 

These  two  apsects  are,  of  course,  related.  For  example,  if  the 
spatial  grid  is  sufficiently  small  to  represent  charge  and  current  densities 
correctly,  then  the  time  step  (Courant  limited  by  the  spatial  grid)  would  be 
small  enough  to  resolve  the  particle  orbit.  However,  the  computational  expense 
associated  with  such  a  brute-source  approach  is  prohibitive.  Thus,  the  desire 
is  to  treat  the  particles  in  some  special  way  (e.g.,  guiding  center,  subcycling, 
etc.)  to  avoid  the  small  spatial  grid  and  time  step. 

To  demonstrate  the  problem  graphically,  we  present  the  results  of 
single  particle  orbits  calculated  with  the  MAGIC  kinematics  algorithm  using 
different  time  steps,  i.e.,  "c/ 1 0 ,  r/4,  ...  . 


2.2  CONVENTIONAL  Kilt  MAT ICS  ALGORITHM. 


The  kinematics  algorithm  in  MAGIC  was  originally  developed  by 
Boris^.  Here,  we  discuss  only  that  aspect  which  relates  to  the  magnetic  field 
term.  Physically,  the  equation  of  motion, 


dp 

dt 


e 

m 


v  x  B  , 


(5) 


t 3 .  P.  Boris,  Relativistic  Plasma  Simulation--Optimization  of  a  Hybrid  Code, 
Proceedings  of  the  Fourth  Conference  on  Numerical  Simulation  of  Plasmas, 

Naval  Research  Laboratory  (1970),  p.  3. 


should  result  in  a  rotation  in  the  plane  transverse  to  B  through  an  angle. 


P  = 


-e<5t  B  , 
ym 


(6) 


without  change  in  the  magnitude  of  the  momentim  vector.  The  procedure  for 
calculating  the  rotation  in  MAGIC  is  illustrated  in  Figure  1.  First,  an 
intermediate  momentun  vector  is  calculated  from  the  vector  addition, 


pb  =  pa  +  flPa  x  0  ’ 

where 


fl 


tan 


-e<5t  B 
2ym 


in  MAGIC,  use  is  made  of  the  expansion, 

tan  p/2  *  p/2  +  0.31755( P/2) 3  ♦  0.2033( P/2) s. 


(7) 


(8) 


(9) 


to  calculate  the  geometrical  factor,  f*.  Finally,  the  rotation  is  completed 
with  the  calculation, 


p=p+f2p.  xB,  (10) 

c  a  b 

where 

,  2f i  (11) 

f2  =  - -  _-2- 

i  +  h  |b| 
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2.3 


CALCULATED  ORBITS. 


From  the  relations  above,  it  is  easy  to  derive  calculated  particie 
trajectories  for  certain  pathological  cases,  e.g.,  6t  =  t,  where  x  is  the  Larmor 
period.  However,  to  demonstrate  this  visually,  we  have  exercised  MAGIC  using 
the  POPULATE  and  PRESET  commands  to  obtain  trajectory  plots  of  calculated  single 
particle  orbits.  The  test  problem  involves  a  1  MeV  electron  in  a  constant 
magnetic  field,  B  =  20T.  The  theoretical  radius  of  gyration  and  Larmor  period 
are 


Ymv  ,  -4 

r  =  =  2.56x10  m 

and  (12) 

x  =  =  5.65x10* 12sec, 


respectively . 

Figure  2  illustrates  particle  trajectory  plots  computed  using  six 
different  values  of  the  time  step,  6t.  Note  that  Figure  2a,  with  ooc6 1  = 
2tc/100,  appears  to  be  a  perfect  circular  orbit.  For  more  useful  time  steps, 
i.e.,  u>c6t  =  2it/10  or  2n/4,  regular  geometric  figures  are  produced  in  which 
the  number  of  sides  is  equal  to  the  time  step  denominator  (see  Figure  2b  and 
c).  For  still  larger  values  of  time  step,  the  orbit  degenerates  into  a  straight 
line;  i.e.,  the  particle  simply  oscillates  back  and  forth.  The  apparent 
oscillation  frequency  is  inversely  proportional  to  the  time  step,  6t,  and  is 
unrelated  to  the  Larmor  period.  Clearly,  this  is  a  poor  representation  of  the 
physical  orbit. 
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(a)  ajc6t  = 


27T/100 


(b) 


wc6t  =  2tt/10 


(e)  w  6t  =  2tt  (f)  a)  6t  =  4tt 

V  C 

Figure  2.  Calculated  particle  orbit  vs.  time  step. 


SECTION  3 

SUBCYCLING  ALGORITHM  IN  MAGIC 


3.1  SUBCYCLING  OPTIONS. 


The  basic  idea  behind  particle  subcycling  is  that  the  particle 
kinematics  uses  a  reduced  time  step;  i.e.,  one  smaller  than  the 
electromagnetic  time  step.  Of  course,  there  are  many  possible  variations 
in  a  subcycling  algorithm;  e.g.,  whether  to  use  the  same  kinematical 
forces  throughout  the  subcycle,  whether  to  treat  electrons  separately  from 
ions,  and  so  forth.  To  test  the  effect  of  such  variations,  the  algorithm 
implemented  in  MAGIC  contained  the  options  listed  below.  (The  flag, 
iscopt,  is  an  input  data  flag  set  by  the  user  in  the  KINEMATICS  command.) 


Subcycie  Option 


no  subcycles 

retain  initial  forces  and  final  currents 
subcycie  forces  and  retain  final  currents 
retain  initial  forces  and  subcycle  currents 
subcycle  both  forces  and  currents 


3.2  FLOW  DIAGRAMS. 


An  overview  of  the  particle  processing  section  of  CONTROL  is 
presented  in  Figure  3.  First  a  check  Is  done  for  particle  creation 
surfaces  (NCRE).  If  particle  creation  has  been  requested,  then  a  flag 
(LSFP)  is  checked  to  ascertain  whether  the  particle  groups  should  be 
processed  in  this  time  step.  If  so,  then  a  check  is  made  (NPRS)  for 
whether  analytically  determined  forces  are  to  be  used  to  move  particles. 
If  not,  then  the  average  vaiue  of  the  electromagnetic  fields  at  the  fuil 
grid  points  is  calculated  by  AFIELOS.  Next,  the  particle  creation  is 
performed  in  CREATOR.  A  test  is  performed  for  the  existence  of  partic'e 
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groups  (NGTI).  Particles  in  MAGIC  are  processed  in  groups  of  192 

particles.  A  call  is  made  to  STORAGE  for  the  next  particle  group  to  be 
processed.  A  loop  over  particle  subcycles  has  been  placed  inside  of  the 
loop  over  particle  groups.  The  particle  kinematics  time  step  is 

determined  by  dividing  the  electromagnetic  time  step  by  the  subcycle 
frequency  (NFSC). 

Figure  A-  is  a  flow  diagram  showing  how  the  subcycle  options  have  been 
implemented  for  particle  forces  and  kinematics.  A  test  on  the  subcycle 
option  flag  (LFSC)  and  the  subcycle  loop  counter  (1TSC)  Is  made  to 

determine  if  particle  forces  are  to  be  subcycled.  When  LFSC=1  or  LFSC=3 
and  ITSC  is  greater  than  1,  then  the  calls  to  the  subroutines,  PREDICT 

which  calculates  the  electromagnetic  forces  at  the  particle  position,  or 
PRSCRIB,  which  determines  analytical  values  for  the  electromagnetic 
forces,  are  calculated  for  the  first  subcycle  and  bypassed  for  subsequent 
subcycies. 

For  LFSC=2  and  IFSC=4,  the  forces  are  calculated  for  every  subcycle. 
The  corrector  flag  (LCOR)  is  initialized  next.  The  subroutines,  K1NEMAT, 
which  computes  particle  trajectories,  and  ROTATOR,  which  transforms 
particle  coordinates  and  momenta,  are  called  for  every  subcycle.  If  the 
predictor  flag  has  been  set  (LPRE=1),  then  the  corrector  flag  is  set 
(LC0R=1)  and  CORRECT,  which  corrects  the  electromagnetic  forces  for  the 
new  particle  positions  through  a  modified  Euler  method,  or  PRSCRIB,  which 
computes  analytical  forces,  is  called  and  K1NEMAT  and  ROTATOR  are  called 
again.  After  the  particle  trajectories  have  been  computed,  then,  if 
requested,  the  trajectories  are  written  to  a  file  for  plotting  by  TRAOECT. 

Figure  5  is  a  flow  diagram  showing  how  the  subcycle  options  have  been 
implemented  for  current  density  calculations.  When  the  particle 
kinematics  are  completed,  the  subcycle  option  flag  is  evaluated.  If 
LFSC=1  or  LFSC=2,  then  SUBCYCl.  is  called,  which  for  these  options  saves 
the  initial  particle  coordinates  and  momenta  of  the  first  subcycle  for 
calculation  of  the  current  density  on  the  last  subcycle.  The  weighting 
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Figure  5.  Current  density,  charge  density,  subcycle, 
and  symmetry  algorithms. 
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fractions  for  current  density  for  the  initial  coordinates  and  momenta  are 
recalculated  in  SUBCYCL  for  LFSC=2,  because  the  subcycling  of  particie 
forces  has  overwritten  the  information.  The  initial  particle  data  is 
replaced  by  the  final  particle  data,  and  the  program  control  goes  to  the 
end  of  the  subcycle  loop  until  the  last  subcycle.  If  LFSC=3  or  LFSC=4, 
the  current  density  is  calculated  for  every  subcycle.  First,  any  symmetry 
boundaries  are  accounted  for  in  SYMMETR,  and  then  the  current  density  is 
calculated  in  3F1EL0S.  A  call  to  SUBCYCL  is  made  after  3F 1ELDS  to  swap 
final  particle  data  tor  initial  particle  data  and  to  calculate  current 
density  weighting  fractions  for  LFSC=3,  since  particie  forces  are  not 
computed  every  subcycle.  CONTROL  then  branches  to  the  end  of  the  loop 
until  the  last  subcycle,  where  particle  data  is  written  to  disk  for  phase 
space  and  movie  plotting  by  subroutines  PHASPAC  and  MOVIE,  If  requested. 
Killed  particles  are  accounted  for  In  SURVIVE  and  charge  density  at  full 
grid  points  is  computed  in  DENSITY. 

When  the  subcycle  loop  and  the  particle  group  are  completed,  then 
J80UNDS  takes  care  of  the  effects  of  boundary  conditions  on  current 
density. 

3.3  KINEMATICS  COMMAND. 

Function:  Specifies  the  particle  kinematics  algorithm. 

Format: 

KINEMATICS  itmult  lrel  apcopt  a3dopt  iscopt  nfsc  / 


Arguments: 

Itmult  -  Particle-to-field  time  step  ratio  (integer). 

lrel  -  Relativistic  flag  (Integer) 

apcopt  -  Predictor-corrector  option  (alpha). 

=  YES  or  NO. 

a3dopt  -  Three-dimensional  option  (alpha). 

=  YES  or  NO. 
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iscopt  -  Particle  subcycle  option  (integer). 

=  0,  no  subcycle  (default). 

=  1,  retain  Initial  forces  and  final  currents. 

=  2,  subcycle  forces  and  retain  final  currents. 

=  3,  retain  initial  forces  and  subcycle  currents. 

=  4,  subcycle  both  forces  and  currents, 
nfsc  -  Subcycle  frequency  (integer). 

=  1,  (default). 

Description: 

The  KINEMATICS  command  is  used  to  specify  the  particle 
algorithm.  The  particle-to-f ield  time  step  ratio  determines  whether 

particle  kinematic  calculations  are  performed  for  every  field  calculation 
or  some  multiple  of  the  field  calculation.  Only  relativistic  kinematics 
are  presently  available  in  MAGIC:  the  relativistic  flag  input  is  simply 
ignored  by  the  code.  The  predictor-corrector  option  can  be  used  to 
estimate  forces  on  a  particle  throughout  the  duration  of  the  time  step. 
It  recalculates  the  particle  trajectory  based  on  that  correction.  This 
requires  two  cycles  through  the  kinematics  algorithm  and  is  therefore  more 
expensive.  The  full  three-dimensional  kinematics  option  is  usually 

selected,  but  the  simpler  two-dimensionai  kinematics  can  be  used  when  the 
TE  mode  is  not  present  or  when  motion  in  the  third  dimension  is  not 
relevant. 

The  particle  subcycle  option  enables  the  particle  kinematics 
calculation  to  be  performed  in  a  multiple  of  steps  or  subcycles  for  each 
calculation  of  the  fields.  The  options  Include  no  subcycling  of  the 
kinematics  (which  is  the  default),  subcycling  only  the  particle  trajectory 
calculations,  subcycling  the  weighted  forces  at  particle  position  along 
with  the  particle  trajectories,  subcycling  the  current  density  calculation 
along  with  the  trajectories,  subcycling  the  forces  and  the  current  density 
along  with  the  particle  forces. 


14 


SECTION  4 

SUBCYCLING  CYCLOTRON  INSTABILITY 


4.1  OVERVIEW. 

This  section  analyzes  an  instability  we  have  termed  the  Subcyciing 
Cyclotron  (SC)  instability.  The  analysis  Is  restricted  in  its  validity  to 
systems  in  which 

co  /  co  >1  .  (13) 

c  p 

where  up  and  He  are  the  electron  plasma  and  cyclotron  frequencies, 
respectively.  It  can  be  shown,  however,  that  this  condition  will  occur  in  many 
systems  for  which  subcycling  would  otherwise  be  useful.  The  SC  instability  is 
therefore  a  serious  problem  and,  if  methods  cannot  be  found  to  eliminate  it, 

this  may  seriously  limit  the  utility  of  subcycling  in  simulations . 

A  derivation  of  the  SC  instability  will  be  done  using  the  fluid  theory 

of  cyclotron  waves.  It  will  be  shown  that  for  the  conditions  of  Equation  (13) 

and  the  limit  of  exact  equations  of  motion  and  fields,  that  only  stable  electron 
cyclotron  waves  will  result.  The  electric  fields  due  to  charge  fluctuations 
will  be  only  a  perturbation  on  the  electron  motion,  and  the  wave  frequency  will 
be  nearly  the  cyclotron  frequency. 

In  the  case  of  subcycling,  however,  the  perturbing  electric  fieid  will 
no  longer  lead  to  a  stable  wave.  The  effect  of  subcycling  will  lead  to  an 

effective  "phase  lag"  of  the  electrostatic  perturbation  on  the  electron  orbit. 
This  phase  lag  will  lead  to  a  growing  amplitude  of  the  electric  fields  and  the 
electron  energy.  This  instability  will  be  fastest  growing  for  the  condition, 

a)  6t  s  it  ,  (14) 

c 

where  6t  is  the  electromagnetic  time  step.  The  SC  instability  will  be  clearly 
distinguishable  from  a  physical  instability  In  that  it  will  lead  to  a  growth  in 
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total  system  energy.  A  physical  instability,  conversely,  will  conserve  total 
energy,  merely  transforming  it  from  one  degree  of  freedom  to  another. 
Therefore,  the  instability  will  then  lead  to  an  unphysical  heating  of  electrons 
and  a  deterioration  of  energy  balance  in  simulations  employing  subcyciing. 

The  results  of  simulations  confirming  the  basic  features  of  the 
analytic  theory  will  be  shown.  These  simulations  were  performed  for  charge 
neutral  systems  but  are  fully  relevant  to  the  non-neutral  case.  The  simulations 
show  a  qrowing  instability  in  the  system  associated  with  the  degree  of 
subcyciing. 

4.2  FLUID  TICORY  OF  CYCLOTRON  WAVES. 

We  will  use  the  fluid  description  of  the  electrons  found  in  a  text  by 
Chen^.  Following  this  treatment,  we  can  derive  the  equations  governing 
electron  cyclotron  waves  In  a  cold  (monoenergetic)  electron  plasma  with  a 
neutralizing  ion  background. 

We  begin  with  the  linearized  equations  for  motion,  continuity,  and 
electric  field,  respectively.  These  are 

md^V  i  =  -e(V  i\8(j+E  A )  , 

a  n;  =  -n„Div(V.)  ,  (15) 

t 


and  Ei  =  -4n  enx  , 

where  Ei ,  n;,  and  Vi  are  the  first  order  perturbed  electric  field,  electron 
density,  and  electron  velocity  respectively,  and  V0,  BQ,  and  nQ  are  the 
steady-state  velocity,  magnetic  field,  and  electron  density,  respectively.  We 
assumed  perturbations  of  the  form,  exp  (kx-wt).  We  break  the  equation  of  motion 
up  into  its  components, 


^Francis  F.  Chen,  Introduction  to  Plasma  Physics,  Plenum  Press  (1974),  p  90. 
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- ioumV x  =  -eE  -  eVyB0  , 

-i<4nVy  =  -eVx8o  ,  (16) 

and  -itijmVz  =  0 

We  solve  this  system  for  Vx  to  obtain 

-iU)mVx  =  e£x  +  eB0  if£o  Vx  ,  (17) 

which  becomes 


oj2  eE 


oj  moo 


where  wc  is  the  electron  cyclotron  frequency.  We  can  now  solve  for  Ex  in 
terms  of  Vx  to  close  the  system.  We  have  from  the  continuity  equation, 


-iojn.  =  innkV 

x 

which  becomes 


"1 


nokVx 

w 


We  then  use  Poisson's  equation  to  write 


(19) 


(20) 


Ex 


-4-rtenj 


ik 


which  then  becomes 


E 


x 


-4nen0Vv 

—  i——  ..A 

iu) 


(21) 


(22) 
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This  expression  is  inserted  back  into  Equation  (18)  to  yield 


co 


x 


CO 


2 


V 

x, 


(23) 


2  2 

where  Wp  =  4ne  n0/m  and  cop  is  the  electron  plasma  frequency.  This  gives 
a  dispersion  relation  for  cold  electron  cyclotron  waves, 


2  2  2 

co*  =  co  z  +  co  *  .  (24) 

c  p 

The  frequency  for  these  waves  is  real  and  therefore  the  waves  are  stable.  For 
2  2 

the  case  <op  / coc«1 ,  we  find  that  the  waves  occur  at  nearly  the  cyclotron 
f requency , 

2 

1  U 

U>  »  w  fw  I  -£  )  .  (25) 

c  2  co2 

c 

%.3  INSTABILITY  DUE  TO  SUBCYCLING. 

We  have  shown  that  the  cyclotron  motion  of  electrons  in  a  plasma  could 

lead  to  stable  electrostatic  oscillations  and  that,  in  the  limit 

2  2 

Wp  /Wc  ^  ’  these  waves  have  a  frequency  nearly  that  of  the  electron 
cyclotron  frequency.  We  now  consider  what  will  happen  when  the  particle 
equations  of  motion  are  advanced  with  a  time  step  much  shorter  than  the 
electromagnetic  time  step.  Defining  a  parameter, 

Se  =  u>c6t  ,  (26) 


we  will  show  that  subcycling  will  perturb  the  cyclotron  waves  and  cause  them  to 
become  unstable,  with  a  growth  rate  depending  on  Se.  The  instability  will 
occur  because  the  electrostatic  field  develops  an  effective  "phase  lag"  relative 
to  the  particle  velocity. 
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The  phase  lag  occurs  because  the  electromagnetic  time  step  defines  a 
period  over  which  the  electric  field  due  to  charge  fluctuations  will  be  frozen 
in  time.  Thus  an  electron  experiences  fields  due  to  the  charge  distributions 
existing  at  the  end  of  the  last  time  step  rather  than  the  instantaneous  charge 
distribution.  For  a  sinusoidal  electric  field,  this  can  be  approximated,  for  a 
small  enough  time  step,  as  a  phase  lag.  This  is  shown  graphically  in  Figure 
6.  Therefore,  for  (XSe<2n,  the  onset  of  subcycling,  we  can  approximate  the 
phase  lag  as  a  phase  factor, 


P  = 
e 


iS  /  2 
e  e 


(27) 


This  approximation  is  based  on  the  relationships  illustrated  in  Figure  6.  The 
phase  factor  can  be  multiplied  times  the  electric  field  from  Equation  18,  to 
yield 


V 

x 


ieE 

mw 


iS  U 
e  e 


(28) 


The  electric  field  is  found  from  Equation  22.  However,  because  the  electron 
velocity  itself  is  advanced  in  discretized  steps,  the  electron  magnetic  force 
will  experience  a  phase  lag  of 


P 

m 


iS  /2 
e  m 


(29) 


where  Pm  =  wc6tp  with  6tp  being  the  time  step  for  advancing  particle 
motion.  This  phase  lag  for  magnetic  forces  occurs  because  the  Lorentz  force  on 
the  electrons  depends  on  their  velocity.  The  velocity  used  to  calculate 
magnetic  forces  is  (in  the  absence  of  predictor-corrector  algorithms)  the 
velocity  found  at  the  end  of  the  previous  time  step.  From  Section  2,  we  have 
demonstrated  useful  values  of  this  phase  lag  to  lie  in  the  range, 

(XS  Oi/4  .  (30) 
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6t 


(a)  Approximation  of  physical  field  over  interval  6t 
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(b)  Approximate  sinusoidal  field 


(c)  Physical  and  approximate  sinusoidal  fields 


Figure  6.  Electric  fields  vs.  time. 
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The  case  of  most  interest,  u)p  /uic  «1 ,  yields 

u)  =  a)  {l  +  —  — £  [cos(Se/2  -  it/8)  +  i  sin  (S J2  -  n/8)J}  .  (37) 

c 

The  imaginary  part  of  the  frequency  is  the  growth  rate  and  has  the  value 

Im  (go)  =  '  _P  sin  (S  /2  -  n/8)  .  (38) 

2  OJ  e 

c 

This  formula  for  the  SC  instability  growth  rate  is  an  approximation  valid  only 
for  0<Se<2n.  However,  we  can  test  the  scaling  by  simulation.  The  scaling 

predicts  stability  for  Se=u/4  and  instability  for  Se>m/4.  A  peak  of 
instability  should  be  reached  at  Se  =  5u/4  with  a  decrease  thereafter. 

4.4  SIMULATION  RESULTS. 

To  test  the  theoretical  model  of  the  SC  instability,  several 

simulations  were  run  with  differing  values  of  Sg.  The  simulations  consisted 

of  a  charge  neutral  plasma  of  electrons  and  ions  with  a  number  density  of  n  = 
1 U  3 

10  /m  in  a  magnetic  field  of  0.5  Tesla.  The  electrons  were  all  given  the 

Q 

initial  random  speed  of  2.7  x  10  m/sec.  This  gives  a  ratio  of  plasma  frequency 
to  cyclotron  frequency,  ojp  /u>c  =  10*  .  The  problem  was  run  for  a  time  of 

10* 10  seconds  or  roughly  100  cyclotron  periods. 

The  distinguishing  feature  of  the  instability,  which  differentiates  it 
from  a  physical  instability  at  the  same  frequency,  is  the  fact  that  energy  is 
not  conserved.  In  a  physical  instability,  particles  give  up  energy  to  a  wave 
field  so  that,  energy  changes  form  but  Is  globally  conserved.  In  the  SC 

instability,  both  particles  and  wave  fields  gain  energy  at  the  same  time,  so 

that  energy  is  not  globally  conserved  but  instead  grows. 

In  Figure  7,  the  instability  is  seen  in  terms  of  global  energy 

growth.  The  total  system  energy  In  fields  ar.d  particles  is  plotted  versus 
time.  The  stable  case,  Se  =  rc/4,  shows  a  flat,  straight  curve  of  energy 
versus  time,  indicating  no  energy  growth.  In  the  very  unstable  case,  Sg  = 

it/2,  Figure  7  shows  that  the  total  system  energy  grows  rapidly,  violating 
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energy  conservation.  The  slopes  of  the  energy  versus  time  curves  in  Figure  7 
are  proportional  to  the  growth  rates  of  the  SC  instability.  The  energy  growth 
rates  for  early  time,  when  linear  theory  is  valid,  are  shown  in  Figure  8.  As 
can  be  seen,  the  peak  growth  rate  is  in  the  neighborhood  of  Se  =  it. 

The  behavior  at  late  time  shown  in  Figure  7  indicates  that  some 

nonlinear  effects  occur  to  slow  down  the  growth  rate.  These  effects  occur  when 

the  wave  fields  are  no  longer  a  small  pertubation  on  the  particle  motion  but 

begin  to  dominate  it.  Trapping  of  particles  in  the  wave  fields  is  a  physical 

effect  which  can  occur  to  slow  the  growth  of  waves. 

4.5  DISCUSSION. 


It  has  been  shown  theoretically  and  demonstrated  with  simulations  that 
subcycling  will  lead  to  an  instability  for  the  case,  Up  /wc  Cl  .  This 
instability  will  lead  to  unphysical  growth  of  both  particle  and  wave  energies. 
The  loss  of  energy  conservation  in  simulations  is  then  one  sign  of  SC 
instability.  The  regime  in  which  the  SC  instability  was  demonstrated  is  a 
restricted  one,  corresponding  to  finite  densities  and  weak  subcycling.  However, 
as  has  been  discussed  in  Section  1,  subcycling  is  only  useful  in  cases  where 
wc6t>1.  For  the  problem  of  interest,  a  high  current  diode  with  regions  of 
magnetic  insulated  sheaths,  we  will  have  the  condition 

Up  =  u)c  •  (39) 

However,  in  those  regions,  subcycling  cannot  be  implemented  without  loss  of 
accurate  particle  dynamics.  This  is  due  to  the  requirement, 

wp6t  <  1,  (40) 


which  must  hold  if  the  dynamics  of  particles  in  electrostatic  fields  is  to  be 
simulated  accurately.  If  o>p»u)c,  then  subcycling,  which  requires  u>c6t>1, 

cannot  be  in  effect.  The  only  regime  of  accurate  electrostatic  particle 
dynamics  and  useful  subcycling  is  therefore 


oo  6t 
c 

<o  5t 
P 


U) 

c 

u 

p 


>  1, 


(41) 
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which  is  the  region  of  SC  instability.  Such  regions  will  occur  outside  the 
sheath  regions  of  the  problem.  Therefore,  at  least  for  magnetically  insulated 
systems,  the  only  regions  where  subcyciing  would  be  useful  are  those  for  which 
it  is  unstable. 
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